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1.0   Introduction 

The  problem  of  constructing  interpolation  functions  for  sets  of 
scattered  data  in  two  independent  variables  has  been  treated  in  many  papers. 
The  survey  paper  on  approximations  to  multivariate  functions  by  Schumaker  [20] 
contains  extensive  references.   Other  survey  papers  are  by  Barnhill  [3] 
and  Sabin  [19].   The  present  author  has  surveyed  and  tested  a  number  of 
algorithms  for  solution  of  the  problem  [10] ,  [11] .   Conceptually  the 
problem  is  quite  simple:   Given  points  (x,  ,y,  ,  f ,  ),  k  =  1,  ...,  N,  with 
distinct  (x,  ,y,  ),  construct  a  function  F(x,y)  so  that  F(x,  ,y  )  =  f.  , 
k  =  1,  ...,  N.   Generally  one  wants  to  have  a  smooth  interpolant,  F(x,y), 
in  the  sense  that  low  order  partial  derivatives  are  everywhere  continuous. 
This  is  complicated  for  large  sets  of  data  by  the  fact  that  the  interpolant 
(in  a  practical  sense,  to  be  computable)  must  be  local,  so  that  its  value 
at  some  point  (x,y)  depends  only  on  (x,  ,y,  ,f,  )  values  for  which  (x,  ,y.)  is 
"close"  to  (x,y).   A  general  framework  for  a  class  of  such  methods  is 
given  in  [9],  and  we  will  discuss  it  briefly  in  Section  2. 

While  a  large  number  of  ideas  have  been  proposed  for  solution  of 
the  problem,  a  much  smaller  number  of  working  computer  programs  are  readily 
available.   These  include:   (1)  a  method  based  on  finite  element  functions 
defined  over  triangles,  due  to  Akima  [1],  [2],  a  version  of  which  is 
available  in  the  IMSL  library  under  the  name  IQHSCV,  (2)  a  program  based 
on  a  similar  idea,  due  to  Lawson  [14],  (3)  two  programs  based  on  weighted 
local  approximations  by  quadratic  functions,  due  to  Franke  and  Nielson  [12]. 
A  program  by  Little  [15],  and  another  by  Nielson  [18],  both  based  on  finite 
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element  functions  defined  over  triangles,  will  probably  be  available 
by  press  time.   All  of  these  programs  have  been  tested  by  the  author 
[10],  [11],  and  most  perform  adequately  in  a  variety  of  cases;  None  of 
them  seems  to  have  a  clear  edge  over  all  the  others,  or  to  be  entirely 
satisfactory.   For  certain  applications,  each  has  its  good  points.   The 
choice  of  a  method  for  most  users  will  be  based  on  subjective  criteria 
which  vary  from  person  to  person  and  application  to  application.   It  is 
not  surprising  this  is  the  case  in  two  variables  since  it  is  also  the  case 
in  one  variable,  although  perhaps  to  a  lesser  extent. 

The  purpose  of  this  paper  is  to  document  an  alternative  scheme  which 
performed  comparably  well  in  the  previously  mentioned  tests.   In  this  way 
the  computer  program  will  be  brought  to  the  attention  of  potential  users, 
who  may  request  it  from  the  author.   The  method  will  supplement  the  current- 
ly available  codes  in  that  it  is  based  on  a  different  approach.   It  is 
anticipated  that  it  will  find  application  and  approval  in  a  variety  of  areas. 

The  theoretical  background  for  the  method  is  discussed  in  Section  2, 
while  details  of  the  program  are  outlined  in  Section  3.   Section  4  gives 
information  concerning  usage  of  the  program.   Several  examples  are  given 
in  Section  5,  including  perspective  plots  of  some  surfaces  generated  by  the 
program.   Included  is  an  example  of  how  the  variation  of  a  parameter  in 
the  method  affects  the  surface.   General  guidelines  for  choice  of  this 
parameter  are  given,  even  though  the  suggested  value  usually  leads  to 
satisfactory  results.   A  general  discussion  of  the  features  of  the  method 
is  given,  along  with  general  guidance  for  use  of  the  program. 

2.0  Theoretical  Background 

The  general  idea  encompassing  this  scheme  and  many  others  are  given 
in  [9].   Consider  that  the  points  (x^y^f^,  k  =  1,  ...,  N  are  given. 
Briefly,  local  approximations  to  the  data  are  constructed,  and  these  are 


then  blended  together  by  using  weighting  functions  which  form  a  partition 

2 
of  unity  on  R~ .   We  pause  to  give  the  necessary  notation  and  definitions. 

A  set  of  functions,  w.(x,y),  i  ■  1,  ,..,  a  is  said  to  form  a  partition 

m 
of  unity  if  each  w.(x,y)  s  0  and   £  w-(x,y)  =  1.   The  w.  will  be  called 

1  i=l   X  X 

weight  functions.   Let  the  support  of  w.  be  denoted  by  S.  =  closure  {(x,y): 

w.(x,y)  5*  0}.   Let  I.  =  {k:   (x,  ,y,)  £  Interior  (S.)}.   Now  suppose  that  the 

functions  Q.(x,y),  i  =  1,  . . . ,  m,  are  defined  on  S.  and  have  the  property 

that  they  interpolate  the  data  whose  (x,y)  coordinates  are  in  Interior  (S.), 

i.e.,  if  k  £  I.,  then  Q.(x,  ,y,)  =  f .  .   These  functions  Q.  will  be  called 

m 
local  approximations.   We  then  consider  the  function  F(x,y)  =  £  w. (x,y)Q. (x,y) 

i=l   X      X 
Its  properties  are  summarized  in  the  following. 

m 
Proposition.    The  function  F(x,y)  =   £  w. (x,y)Q. (x,y)  has  the  following 

i=l  X      1 
properties: 

(1)  Interpolation;  F(x,  ,y  )  =  f,  ,  k  =  1,  ...,  N. 

(2)  Smoothness;  F(x,y)  is  at  least  as  smooth  as  the  w.  and  Q.,  e.g., 
if  all  of  the  functions  w. ,  Q.,  i  =  1,  ...,  m  have  continuous 
first  derivatives,  so  does  F(x,y) . 

(3)  Local  dependence  on  the  data;  Let  (x,y)  be  fixed,  and  let 

J    =  {i:   w„(x,y)  ?*  0},  then  F(x,y)  depends  only  on  the  (x,  ,y,  ,f,) 

X,  y  X  K   K   K 

points  for  which  k  e(  u   I.)  u(i:   some  0.,  i  £  J    depends  on 

UJ    x  J   J    x,y   ^ 

x,y 

(x.,y.,f.)}.   In  particular,  we  have  F(x,y)  =     I       W. (x,y)Q. (x,y) . 

ieJ 

x,y 

These  properties  are  essentially  observations,  but  form  the  basis  for 
construction  of  appropriate  weight  functions  which  will  allow  easy  determin- 
ation of  the  set  J   .   Our  construction  yields  a  set  of  at  most  four  nonzero 
x,y 

terms  in  the  sum  defining  F(x,y).   This  provides  a  considerably  faster  process 
during  the  evaluation  of  the  interpolant  than  was  possible  in  the  choices 
previously  considered  [9]. 


It  has  been  implied,  but  is  not  necessary  for  the  proposition  to 
hold,  that  many  of  the  weight  functions,  w.  ,  have  finite  support.   In  order 
for  the  local  approximations  Q.(x,y)  to  actually  be  local,  this  will 
likely  be  the  case.   Therefore  we  think  of  weight  functions  whose  support, 

S.,  contains  relatively  few  (x,  ,y,)  points,  and  whose  support  is  often 

2 
local.   In  order  for  the  interpolant  to  be  defined  on  all  of  R  ,  some  sets 

S.  must  not  be  finite,  and  typically  would  contain  points  (x^,y,  )  on  or 

near  the  boundary  of  the  convex  hull  of  the  set  of  points  {(x,  ,y,)}.   With 

this  framework  and  ideas  in  mind,  we  are  ready  to  discuss  the  specific 

details  of  the  algorithm. 


3.0  Details  of  the  Algorithm 


3.1  Choice  of  Weight  Functions 

While  the  choice  of  weight  functions  was  allowed  to  determine  the 

support  regions  in  the  discussion  of  the  previous  section,  it  is  more 

convenient  to  proceed  from  support  regions  to  weight  functions  in  the 

actual  application.   The  use  of  support  regions  which  are  rectangles  have 

specific  advantages  in  terms  of  controlling  the  number  of  support  regions  in 

which  a  particular  point  (x,y)  lies,  as  well  as  simplifying  determination 

of  those  regions. 

Let  n  and  n  be  given  positive  integers,  and  let  finite  values  of 
x      y 

xQ  <  X]_  <  x2  <  ...  <  xn   <  xn  +1  and  yQ  <  y  <  ^    <  .  .  .  <  yn   <  yn  +±   be 

xx  y    y 

given.   For  each  i  =  1,  2,  ...,  n  and  j  =  1,  2,  ...,  n  let  r.  .  = 

x  y  i  j  j 

Cxi-r  xi+i]  *  frj-r  ^j+i]- 

2     3 
Let  H^Cs)  =  1  -  3s  +  2s  ,  the  Hermite  cubic  satisfying  H.,(0)  =  1, 

H^Cl)  =  H'(0)  =  H'(l)  =0.   We  then  define  piecewise  cubics  with  continuous 


first  derivatives,  which  are  nonzero  only  on  two  adjacent  intervals,  and 
satisfy 


vi(x  )  =  5    ,    i,j  =  1,  2,  . ..,  nx 


u.(y.)  =5..     i,j  =1,  2,  ...,  n  , 
J  1    ji  y 


In  particular  we  take 


r  i 


v1(x) 


=  <  H 


X  -  X, 


x„  -  X, 


X  <  X, 


X   <  X  <  X- 


X  s  x. 


r 


v.(x)   =  { 

1 


X  -  X. 

1 


,     X  <  X 


i-1 


1  -  v ,  . (x)    ,    X.  .  <  X  <  X. 

1-1  1-1        1 


x.  ,.  -  X. 

1+1    1 


V 


X.  <  X  <  X. , , 

1       1+1 


X  >  X 


i+1 


for  i  =  2,  . . . ,  n  -  1 

X 


,     X  <  X 


vn  (x)  =  • 

X 


n  -1 
x 


1  -  v    ,  (x)   ,    X    t  <  X  <  X 

n  -1  n  -1        n 

X  XX 


,     X  >  X 


The  u.(y)  are  dual.   Then  if  we  define 


W   (x,y)  =  v  (x)u  (y) 

i)j         i    i 


i  =  1,  . . . ,  nx  ,   j  ■  1,  . ..,  n  , 


it  is  easily  observed  that  \J.        has  support  r.  .,  except  for  when  i  =  1  or  n  , 

i>3  i>3  x 

or  j  =  1  or  n  ,  when  the  support  extends  to  OT  in  one  or  both  variables.   We 

denote  the  support  of  W.  .  by  R.  ..  Further,  we  note  that  the  W.    form  a 

1,3      i,J  i»3 

partitxon  of  unity  on  R  . 

Since  the  x.  and  y.  values  give  rise  to  a  grid  on  the  plane,  we  call 
them  grid  values.   The  choice  of  grid  values  x.  and  y.  depend  on  the  data. 
They  may  be  specified  by  the  user,  but  the  default  option  is  for  them  to  be 
determined  by  the  program.   In  the  default  mode  the  user  specifies  a  parameter, 
NPPR  (for  number  of  points  per  region),  the  suggested  value  being  about  10. 
The  program  will  then  determine  the  grid  values  so  that  the  anticipated 
number  of  (x,  ,y.)  points  in  each  rectangle  R.  .  will  be  approximately  NPPR. 

K,   K  X,J 

For  data  which  is  not  somewhat  uniformly  distributed  the  actual  numbers  may  vary 

considerably.   However,  for  most  situations  we  have  encountered,  the  process 

is  quite  adequate. 

An  equal  number  of  grid  values  is  taken  in  each  direction,  i.e.,  n  =  n  . 

x    y 

Because  we  want  NPPR  points  per  rectangle,  each  subrectangle  (x.,x  ..)  x  (y ,  ,y .  A 

1  2  1 

should  have  7-  NPPR  points.   Since  n  =  n  ,  we  want  (n  +1)  —  NPPR  =  N, 
4  x    y  x       4 

1/2 
leading  us  to  take  n   to  be  the  nearest  integer  to  (4N/NPPR)    -  1. 

The  grid  values  x. ,  are  now  chosen  so  that  there  are  approximately  equal 

numbers  of  points  from  the  values  x,  ,  k  =  1,  ...,  N  in  each  interval 

r~  A 

(x  ,  x   ;.   Specifically,  let  x,  denote  the  x,  values  arranged  in  non- 
decreasing  order.   Consider  the  points  (0,x..),  (l,x9),  ...,  (N-l,  O  ,  and 
let  g(t)  be  the  piecewise  linear  interpolant  for  them.   Divide  the  interval 

(0,N-1)  into  n  +  1  subintervals  of  length  A  =  — ~ -   .   The  values  of  x.  are 
x  n  +1  1 

x 

determined  by  taking  them  to  be  the  values  of  g  at  the  endpoints  of  the 

subintervals,  i.e.,  x  =  g(i  ),  i  =  0,  1,  2,  ...,  n  +  1.   The  y.  are  determined 

x  3 

in  dual  fashion.   This  process  results  in  the  grid  values  and  hence  the 
rectangles,  being  symmetric  if  the  (x^.y.)  points  are  symmetric.   In  addition, 
the  relative  positions  of  grid  values  are  unchanged  by  linear  displacements 

6 


and  stretching  in  each  variable. 

When  chosen  in  the  above  fashion,  the  location  of  the  lines  is  not 
a  local  process  in  the  sense  that  insertion  of  an  additional  point  will 
change  the  boundaries  of  all  of  the  rectangles.   While  one  could  argue  that 
then  the  scheme  is  not  local,  we  take  the  view  that  the  idea  of  local 
determination  of  the  interpolant  is  most  important  in  the  evaluation  phase. 
The  determination  of  parameters  in  the  scheme  (here,  the  x.  and  y.)  may  be 
a  global  process.   Of  course,  if  the  user  specifies  the  grid  values,  he  will 
likely  be  using  a  global  process  to  choose  them. 
3.2   Choice  of  Local  Approximations 

The  only  constraint  on  the  local  approximations  is  that  they  interpolate 
the  appropriate  points,  and  that  they  have  continuous  first  derivatives 
(at  least)  to  assure  a  smooth  interpolant.   In  the  previously  mentioned  tests 
conducted  by  the  author,  a  number  of  global  interpolation  schemes  for  scattered 
data  were  considered.   In  principle,  any  of  these  might  be  used.   The  choice 
here  was  made  for  two  reasons,  (1)  the  method  scored  very  well  in  the  tests, 
and  (2)  the  method  has  an  elegant  and  well  developed  mathematical  theory  which 
also  has  direct  application  to  some  engineering  problems. 

The  local  approximations  used  in  this  algorithm  are  the  thin  plate  splines 
first  mentioned  by  Harder  and  Desmairis  [13]  ,  with  theoretical  developments 
by  Duchon  [4-7]  and  Meinguet  [16],  [17].   It  was  first  developed  as  the 
solution  to  the  problem  of  a  thin  plate  which  is  forced  to  pass  through 
certain  points  (the  interpolation  points)  by  application  of  point  loads.   For 
our  purposes  it  is  sufficient  to  know  that  the  approximation  is  of  the  form 

r     2 
Q(x,y)  =  l     A^dk  log  dfc  +  a  +  bx  +  cy  , 

kel 

2  2 

where  I  =  {k:   Q  is  to  take  on  the  value  f,  at  Cx,  ,y,)},  and  d,  =  (x  -  x,  )   + 

2 
(y  -  vO    •      Ttie  coefficients  A,  ,  and  a,  b,  and  c  are  determined  by  the  linear 


system  of  equations 

2 
I     Akdk  log  dk  +  a  +  bx  +  cy 

kel 


(x,y)  =  (xi,yi) 


=  f.  ,   iel 

1 


kel  K 

I  Vk  =  ° 

kel  ^   k 

kel  k  fc 

The  geometric  effect  of  the  last  three  equations  is  to  suppress  all  terms 
in  the  approximation  which  grow  faster  than  linear  as  distance  from  the 
interpolation  points  is  increased.   A  linear  system  of  order  equal  to  the 
number  of  interpolation  points  plus  three  must  be  solved.   To  be  nonsingular 
there  must  be  at  least  three  noncollinear  points  among  the  (x,  ,y,),  kel. 
The  system  is  symmetric,  but  not  positive  definite.   While  an  equation 
solver  designed  for  such  systems  could  be  used,  we  have  found  a  general 
purpose  solver,  the  DECOMP/SOLVE  subroutines  of  Forsythe,  Malcolm,  and  Moler  [8] 
has  given  more  reliable  results. 

It  is  easily  observed  that  the  local  interpolant  has  continuous  derivatives 
of  all  orders  except  at  the  data  points,  (x,  ,y,  ) ,  where  a  logarithmic 

rC    iC 

singularity  occurs  in  the  second  derivatives.   The  interpolant  has  linear 
precision,  that  is,  if  the  (x,  ,y,  ,f  )  points  all  lie  on  a  linear  function,  the 
interpolant  will  reproduce  it. 

While  the  thin  plate  spline  is  invariant  with  respect  to  scaling,  trans- 
lation, and  rotation  (not  all  of  this  is  obvious),  the  condition  number  of 
the  coefficient  matrix  for  the  system  of  equations  is  dependent  on  scaling. 
To  minimize  difficulties  with  that,  and  to  remove  the  effect  of  scaling  the 

variables  by  different  amounts,  we  transform  each  rectangle  r.  .  onto  the  unit 

2 
square  [0,1]  ,  before  the  local  approximation  Q.  .is  computed. 

i»  J 
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It  remains  to  specify  the  process  for  determining  the  points  (x,  ,y,  ,  f ,  ) 

to  be  interpolated  by  the  thin  plate  spline  local  approximation.   Experience 

has  shown  that  it  is  advantageous  to  include  more  points  than  is  necessary, 

i.e.,  (x,  ,y  )  which  are  outside  of  R.  ..   This  tends  to  yield  a  better  trans- 
it.  iC  1,  J 

ition  between  local  approximations  than  when  only  necessary  points  are  included, 

2 
The  set  of  (x,  ,y,  )  points  transformed  into  the  rectangle  R  =  [-.1125,  1.1125] 

2 
by  the  transformation  taking  r.  .  onto  [0,1]   are  included.   Let 

1  >  J 

\   "  Xi-1       7k  "  yi-l 
I.    =  (k:   -^ M^—  ,  -^ :LzL-    \e   R} 

1  x  -  -  x.       y    -  y. 

i+l    i-l      J+l    j~l 

This  gives  the  basic  set  of  interpolation  points  for  the  local  approximation 

Q.  .  associated  with  the  rectangle  R.  . .   Under  certain  conditions  there  may 

be  fewer  than  the  necessary  three  indices  in  I.  ..   When  this  happens,  the 

set  I.  .is  augmented  by  including  as  interpolation  points  the  necesary 
■L  >  J 

number  of  closest  points  to  the  rectangle  R.  .  (  in  the  I     norm),  after  the 
points  (x,  ,y  )  have  undergone  the  transformation  to  the  unit  square.   The 
minimum  number  of  points  per  rectangle  is  a  variable,  MINPTS .   This  has  been 
set  to  3,  but  may  be  increased  if  it  seems  desirable. 

After  the  interpolation  points  for  each  local  approximation  have  been 
determined,  the  local  thin  plate  splines,  Q.  .,  can  be  determined  by 
calculating  the  coefficients.   This  yields 

Q.  .(x,y)  =   I  A.  ,  .  d/2  log  d/  +  a.  .  +  b.  .  x'  +  c.  .  y' 

^x,J   J  k£J     t,j,k  k       k    i,j    i,j       i,j 

where  the  primes  denote  coordinates  and  distance  after  the  transformation  of 

r.  .to  the  unit  square 

3.3   Properties  of  the  Interpolant 

The  overall  interpolant  is  of  the  form 
x   v 


F(x,y)  'I         Y  W  ,(x,y)Q,  .(x,y), 

i-i  A  L*j       ^ 


and  as  noted  previously,  there  are  at  most  four  nonzero  terms  in  the  sum. 

Which  terms  are  nonzero  is  easily  determined.   If  x  ^  x   ,  set  i'  =  n  ; 

n  x 

x 

Otherwise  let  i'  be  the  smallest  index  so  that  x     >  x.   Determine  j1  in 

3.  +1  J 

dual  fashion  from  y  and  the  y.'s.   Then,  the  four  nonzero  terms  have  indices 
(i\j'),(i'+l,j'),  (i'.j'+D,  and  (i'+l,j'+l).   If  i*  -  0  or  i'  =  n^,    the 
terms  involving  i'  or  i'+l,  respectively,  do  not  appear,  and  similarly  if 
j'  =  0  or  j '  =  n  ,  the  terms  involving  j'  or  j'+l,  resoectively ,  do  not 
appear. 

In  addition  to  the  properties  outlined  in  Section  2,  certain  other 
properties  hold.   The  approximation  is  invariant  under  translation  and 
stretching  (independently  in  each  variable) .   It  has  symmetry  with  respect 
to  planes  parallel  to  coordinate  planes  whenever  the  data  has  that  symmetry. 
The  approximation  is  not  invariant  under  rotations,  however,  since  the 
rectangles  depend  on  the  individual  coordinates  of  the  data  points. 

The  approximation  has  continuous  first  derivatives,  and  jump  dis- 
continuities in  the  second  derivatives  across  grid  lines,  as  well  as 
logarithmic  singularities  in  the  second  derivatives  at  the  data  points. 
Plots  of  the  surfaces  generally  appear  to  be  quite  smooth,  however.   Since 
the  local  approximations  have  linear  precision,  the  overall  approximation 
also  has  linear  precision,  i.e.,  if  the  data  lies  on  a  linear  function,  the 
interpolant  is  a  linear  function. 

4.0  The  Computer  Program 

The  overall  hierarchy  of  the  subroutines  is  given  in  Diagram  1,  which 
shows  the  communication  links  between  them.   No  COMMON  is  used  as  array 
sizes  are  problem  dependent  and  thus  specified  by  the  user.   We  briefly 
discuss  them,  mentioning  important  parameters. 
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USER  PROGRAM 

\[ . 

LOTPS 

vr 

> 

f 

\ 

f 

\t  „. 

\( 

GRID 

LOCLIP 

CLOTPS 

EVLTPS 

\ 

t 

w 

VSORTA 
(IMSL) 

DECOMP 
SOLVE 

Diagram  1. 

Calling  Program 

This  program  is  supplied  by  the  user  and  must  supply  the  ^^'^'^^ 
points,  plus  two  arrays  of  points  xO . ,  and  yO .  for  the  grid  of  points 
(xO.,yO.)  at  which  the  interpolant  is  to  be  evaluated.   In  addition,  the 
user  must  supply  two  workspace  arrays,  IWK  and  WK  in  which  information 
calculated  during  preprocessing  (e.g.,  I    ,  x  ,  y  ,  and  coefficients  for 
the  local  approximations),  is  stored,  and  an  array  FO  for  the  returned 
interpolant  values. 

The  amount  of  storage  required  for  the  arrays  IWK  and  WK  is  not 
known  a  priori.   The  estimated  space  required  is  about  6N  for  IWK  and 
about  7N  for  WK.   Table  1  gives  exact  results  for  several  different  sized 
problems  based  on  random  (x,y)  points.   Oddly  distributed  point  sets  may 
result  in  somewhat  more  storage  being  required.   In  any  case,  the  precise 
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number  of  locations  required  is  returned  to  the  calling  program  from 
Subroutine  LOTPS,  the  only  routine  referenced  by  the  user.   If  an  in- 
sufficient number  are  allowed,  an  error  return  occurs. 

Under  the  usual  option,  the  user  specifies  NPPR  >  0,  the  suggested 
value  being  10.   If  the  user  wishes  to  specify  the  grid  lines,  he  may 
do  so  by  setting  NPPR  =  0,  and  then  giving  grid  line  information  in  the 
arrays  IWK  and  WK,  as  explained  in  the  argument  description  for  Subroutine 

LOTPS.   Typically  one  should  take  xn  =  min  x,  ,  x   ..  =  max  x,  ,  and  the  dual 

x 
in  y.   This  is  not  necessary,  although  all  points  to  be  interpolated  should 

lie  in  [x..,  x   . ,  ]    [yA,  y   .  n  ]  .   To  prevent  different  scaling  (internally) 
u   n  +1      U   n  +1 
x  y 

in  the  two  variables,  a  square  grid  covering  the  (x,  ,y,  )  points  could  be 

specified. 

Subroutine  LOTPS 

This  subroutine  provides  the  interface  between  the  user's  program  and 

the  set  of  routines  implementing  the  method.   Generally,  LOTPS  sets  up 

storage  areas  in  the  arrays  IWK  and  WK,  determines  parameters  required  by 

other  subroutines,  and  calls  other  subroutines  to  (1)  generate  the  grid,  if 

necessary,  (2)  determine  the  interpolation  points  for  the  local  approximations, 

(3)  compute  coefficients  for  the  local  approximations,  (4)  evaluate  the 

interpolant  at  a  grid  of  points. 

Subroutine  GRID 

This  subroutine  selects  values  of  x.  and  y.  in  accordance  with  the 

i      3 

discussion  in  Section  3. 
Subroutine  LOCLIP 

This  subroutine  determines  the  interpolation  points  for  each  local 
interpolant,  in  accordance  with  Section  3. 
Subroutine  CLOTPS 

This  subroutine  generates  the  system  of  equations  for  the  coefficients 
of  the  local  approximations  and  calls  an  equation  solver  to  obtain  them. 
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Internally  the  routine  has  a  maximum  of  30  local  interpolation  points. 
This  can  easily  be  altered  by  changing  two  statements,  as  noted  in  the 
program  listing. 
Subroutine  EVLTPS 

This  subroutine  evaluates  the  interpolant  on  the  grid  of  points  specified 
by  the  user.   Use  of  a  grid,  when  that  is  required,  facilitates  the  process 
of  locating  the  rectangle  in  which  a  particular  evaluation  point  is  located. 
The  xO.  and  yO .  values  should  be  in  increasing  order  for  maximum  efficiency. 
Evaluation  time  for  a  grid  of  points  should  be  nearly  independent  of  N  for 
large  N,  which  is  borne  out  by  the  timing  information  in  Table  2. 

5.0   Examples  and  Observations 

Example  1 .   This  example  shows  a  typical  local  approximation  function  and 
is  for  the  data  in  Table  1.   This  function  is  a  "cardinal"  function  for  the 

first  point,  and  as  such  shews  the  effect  of  a  nonzero  value  at  a  single  point 

2 
on  the  interpolant.   The  plotted  surface,  shown  in  Figure  1,  is  over  [0,1]  . 

Example  2.   This  example  is  given  to  show  a  surface  generated  from  100 

randomly  generated  points  with  the  function  value  being  obtained  from  an 

explicitly  given  function.   The  surface  was  generated  using  NPPR  =  10,  has 

rms  error  of  about  .3%,  and  is  virtually  indistinguishable  from  a  plot  of  the 

parent  surface.   It  is  shown  in  Figure  2. 

Examples  3 ,  4 ,  5 .   These  examples  use  the  same  set  of  60  points  lying 

1   19   2 
in  the  square  (  -  -r-g   ,  rr  )  ,  and  chosen  by  a  pseudorandom  number  generator. 

The  function  is  explicitly  given  by  f(x,y)  =  .1  +    .  ,        ^        ,    and  is 

shown  in  Figure  3.   Figures  4,  5,  and  6  show  the  interpolant  over  the  square 

2 
[0,1]   for  NPPR  values  of  6 ,  10,  and  15. 

From  these  examples  we  see  that  NPPR  =  10  works  well,  not  too  much 

difference  is  observed  when  NPPR  is  increased,  but  NPPR  =  6  gives  a  less 


13 


smooth  appearing  surface.   The  smaller  NPPR  is,  the  more  localized  the 
surface  becomes  (although  NPPR  =  4  will  probably  be  the  least  value 
practicable,  and  some  local  interpolants  will  likely  become  planes  due  to 
the  minimum  of  three  interpolation  points  being  reached) .   In  line  with 
this  comment,  very  smooth  surfaces  with  small  gradients  will  probably  be 
amenable  to  larger  NPPR,  while  surfaces  with  large  gradients  may  be  best 
approximated  by  taking  NPPR  smaller,  thus  localizing  the  behavior. 

Table  2  gives  the  results  of  a  series  of  problems  with  various  numbers 

of  points  and  values  of  NPPR.   The  data  points  were  chosen  by  a  pseudorandom 

1   19   2 
number  generator  in  the  square  [  -  tt:  ,  yjr  1  »  and  the  approximation  was 

2 
evaluated  on  a  33x33  grid  (of  1089  points  total)  on  [0,1]  .   We  observe  that 

increasing  NPPR  increases  execution  time  while  decreasing  the  amount  of 

2 
storage  needed.   Preprocessing  time  should  be  about  proportional  to  N  ,  but 

apparently  there  is  a  strong  linear  component  for  small  N.   Preprocessing 

time  should  also  be  about  proportional  to  NPPR,  although  this  is  not  readily 

apparent.   Evaluation  time  should  be  nearly  independent  of  N  for  large  N,  and 

proportional  to  NPPR,  which  is  approximately  shown. 

The  automatic  grid  selection  process  works  well  when  the  data  is 

fairly  uniform  in  (x,y)  and  lies  nearly  in  a  square  region.   If  the  data  is 

very  irregular,  or  lies  in  an  oblong  rectangular  area,  it  will  probably  be 

useful  to  explore  results  with  a  user  specified  grid  to  obtain  better  coverage 

by  rectangles  which  are  not  too  oblong.   Limited  experience  has  been 

accumulated  in  these  situations. 

6.0   How  to  Obtain  this  Program 

A  copy  of  this  program,  written  in  Fortran,  and  including  a  sample 
driver  program,  can  be  obtained  from  the  author.   To  do  so,  send  a  (short) 
1/2  inch  tape  to  the  author  indicating  the  format  desired. 
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0.35  0.35  0.5 

-0.05  0.25  0.0 

0.10  -0.05  0.0 

0.50  0.05  0.0 

0.00  0.90  0.0 

0.30  0.70  0.0 

0.60  0.50  0.0 

0.90  0.00  0.0 

0.40  1.05  0.0 

0.85  0.80  0.0 

1.05  0.20  0.0 

1.10  1.10  0.0 


Table  1:   Data  for  "Cardinal"  function. 
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NPPR 


6 

60 

5 

6 

100 

7 

6 

500 

18 

6 

1000 

25 

10 

60 

4 

10 

100 

5 

10 

500 

13 

10 

1000 

19 

15 

60 

3 

15 

100 

4 

15 

500 

11 

15 

1000 

15 
Table  2 

Time 

Time 

IWKU 

(Preprocessing) 

(Evaluation) 

273 

4.2 

14.3 

506 

7.5 

17.6 

2893 

70.5 

21.2 

6140 

201.5 

21.3 

243 

5.3 

18.0 

427 

9.9 

26.4 

2649 

73.3 

29.8 

5804 

202.6 

31.7 

199 

6.8 

23.0 

373 

12.4 

32.8 

2591 

91.6 

39.3 

5439 

258.5 

40.6 

NWKU 

334 

619 

3596 

7441 

284 

488 

3014 

6565 

224 

414 

2856 

5920 


Storage/Times  for  Various  Size  Problems. 
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Figure  1:   Cardinal  Function 


Figure  2 :   Saddle  Function 
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Figure  3:   Parent  Function 


Figure  4:   NPPR  =  6 


Figure  5:   NPPR  =  10 


Figure  6:   NPPR  =  15 
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